SUBROUTINE WELFARE_GAIN
!Computes welfare gains in terms of permanent increase in total consumption from switching from Economy A to Economy B
! Economy A: to which add welfare gain (e.g. Samuelson)
! Economy B: Better economy (Baseline)

!INSTRUCTIONS:
!1. To compute policy fcns with (na,nb), need to reset nb2, re-run func_gN first and store txt files for policy fcns in folder welfare
!2. Need to call FUNCALLER_GN.F90 beforehand to upload all parameter values

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

REAL(8),DIMENSION(nb,na)::v_matrixA,v0_matrixA,vaut_matrixA
REAL(8),DIMENSION(nb,na)::b_next_matrixA,q_matrixA,q_matrix_nodefA
REAL(8),DIMENSION(nb,na)::x_matrixA,gn_matrixA,xaut_matrixA,gnaut_matrixA
INTEGER,DIMENSION(nb,na)::default_decisionA

REAL(8),DIMENSION(nb,na)::cn_matrixA,cnaut_matrixA,ct_matrixA,ctaut_matrixA,c_matrixA,caut_matrixA
REAL(8),DIMENSION(nb,na)::tau_matrixA,tauaut_matrixA

REAL(8),DIMENSION(61,na)::v_matrixB,v0_matrixB,vaut_matrixB
REAL(8),DIMENSION(nb,na)::vaut_matrixA_w,vaut_matrixA_w_new,v0_matrixA_w_new,v_matrixA_w,v_matrixA_w_new
REAL(8)::dev_v,dev_v2,ctval,cnval,b_next,bpres,p0,q,mstar,a0,gain_c,criteria
REAL(8)::hres,cres,gnres,taures,utilc,utilg,u_value,a_initial,b_initial,aux0

REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::exp_vA,exp_vB,a00,b00 
INTEGER,PARAMETER::nbA=nb2,nbB=nb2
REAL(8)::distAaccess(nbA,na),distBaccess(nbB,na)
REAL(8),DIMENSION(na)::distAaut,distBaut

REAL(8)::a_threshold,a_max1,a_min1,a_next,a_next1,acum2,acum1,acum3,cdf,cdf1,cdf_threshold,scalar
REAL(8)::exp_v_next,exp_v_next_aut,a_fun,value_next,ucost,udef,DNORDF,objective_fun,defval
REAL(8)::c0,gn0,slope,DNORIN
EXTERNAL::a_fun,udef

INTEGER ILEFT, IRIGHT,NINTV,choice_cond,i_a_global,i_b_global,iter,iter2,nbtemp
EXTERNAL DNORIN
INTEGER::i,i_a,i_b,ii_b,iii_b,i_d,i0(1),i_a0,i_b0,i_b0A,i_b0B,i_default_today,i_defaultA,convergence,convergence2
REAL(8)::DCSVAL,gain_c_new,dev_gain ,valuex,dev_v2_old
EXTERNAL::DCSVAL

REAL(8)::v0_fun,gain_c_state(nb),v0_fun_welfare,weight_gain,bgridB(61)
INTEGER::nbaux,index_opt,nd2

EXTERNAL::v0_fun,v0_fun_welfare
REAL(8)::new_value,hval,c01,c02,aahat,aa1,aa0,a01,a02 ,slope1,slope2,lag_lim 
REAL(8)::exp_vA_repay,exp_vB_repay ,v22, valueA,valueB, slope_aux(nb)


open (110, FILE='graphs\b_grid.txt')
DO i=1,nb
    READ(110,*) bgrid(i)
ENDDO

CLOSE(110)

open (210, FILE='graphs_baselineX\v.txt')
open (211, FILE='graphs_baselineX\bgrid.txt')
do i_b = 1,61
    READ(211,*) bgridB(i_b)
    do i_a = 1,na
        READ(210, '(F18.10, X, F18.10, X, F18.10)') v_matrixB(i_b, i_a), &
                    v0_matrixB(i_b, i_a), vaut_matrixB(i_b, i_a)
    end do
end do

CLOSE(210)
CLOSE(211)

open (10, FILE='graphs\v.txt')
open (11, FILE='graphs\default.txt')
open (13, FILE='graphs\b_next.txt')
open (114, FILE='output\welfare\gh.txt')
open (115, FILE='output\welfare\ggn.txt')
open (116, FILE='output\welfare\ghaut.txt')
open (117, FILE='output\welfare\ggnaut.txt')
open (118, FILE='output\welfare\gct.txt')
open (119, FILE='output\welfare\gctaut.txt')

do i_a = 1,na
do i_b = 1,nb
    READ(114, *) x_matrixA(i_b, i_a)
    READ(115, *) gn_matrixA(i_b, i_a)                     
    READ(118, *) ct_matrixA(i_b, i_a)

ENDDO
ENDDO



do i_b = 1,nb
do i_a = 1,na
    i_a_global = i_a
    i_b_global = i_b

    READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrixA(i_b, i_a), &
                v0_matrixA(i_b, i_a), vaut_matrixA(i_b, i_a)
    READ(13, '(F15.11)') b_next_matrixA(i_b, i_a)
    if (vaut_matrixA(i_b, i_a) > v0_matrixA(i_b, i_a)) then
        default_decisionA(i_b, i_a) =2
    else
        default_decisionA(i_b, i_a) =1
    end if
end do
end do

do i_a = 1,na
    READ(116, *) valuex
    xaut_matrixA(:, i_a)=valuex
    READ(117, *) valuex
    gnaut_matrixA(:, i_a)=valuex 
    READ(119, *) valuex
    ctaut_matrixA(:, i_a)=valuex 
enddo    

CLOSE(10)
CLOSE(13)
CLOSE(114)
CLOSE(115)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)

PRINT*,' '
PRINT*,'SOLUTION UPLOADED FOR WELFARE ANALYSIS'
PRINT*,' '

!upload value and policy fcns
PRINT*,'Compute CONDITIONAL consumption gain' 
PRINT*,' '

do i_a = 1,na
do i_b = 1,nb
    a0 = DEXP(agrid(i_a))
    a_initial = agrid(i_a)
    b_initial = bgrid(i_b)
    
    hval = x_matrixA(i_b, i_a)
    cnval = hval**alpha-gn_matrixA(i_b, i_a)
    cn_matrixA(i_b, i_a)= cnval
    
    b_next = b_next_matrixA(i_b, i_a)
    ctval = ct_matrixA(i_b, i_a) 

    IF (mu .NE. 0d0) THEN
        c_matrixA(i_b, i_a) = (omegac*ct_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cn_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        c_matrixA(i_b, i_a) = (ct_matrixA(i_b, i_a)**(omegac))*(cn_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
    
    cnval = xaut_matrixA(i_b, i_a)**alpha-gnaut_matrixA(i_b, i_a)
    cnaut_matrixA(i_b, i_a)= cnval
    ctval = ctaut_matrixA(i_b, i_a)

    IF (mu .NE. 0d0) THEN
        caut_matrixA(i_b, i_a) = (omegac*ctaut_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cnaut_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        caut_matrixA(i_b, i_a) = (ctaut_matrixA(i_b, i_a)**(omegac))*(cnaut_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
ENDDO
ENDDO

OPEN(1111,FILE='output//gain_c_state_min2.txt')

CLOSE(111)

!To apply a_fun, need original value fcns
!Default thresholds computed using model A
vaut_matrix  = vaut_matrixA
v0_matrix    = v0_matrixA
v_matrix     = v_matrixA

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
    ii_b=index_bfeas4(i_a)+1
    nbaux =nb-ii_b+1
    FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
ENDIF         
ENDDO

criteria=1d-5
iter=0
gain_c_state = 0.01d0

a00 = mua/(1d0-rhoa)       !start from unconditional mean
i_a0 = LOCATE2(agrid,a00)

i0 = LOCATE2(bgrid,0d0)
i_b0A = i0(1)
i0 = LOCATE2(bgridB,0d0)
i_b0B = i0(1)

DO ii_b=i_b0A,i_b0A
    i_b0 = ii_b

    IF (MIN(v_matrixB(i_b0B, i_a0),v_matrixA(i_b0A, i_a0)).LE.-500d0) THEN
        gain_c_state(i_b0) = -1d0
        CYCLE
    ENDIF

    convergence = -1
    convergence2= -1

    gain_c = 0.01d0
    dev_gain  = 1d0

    !v_matrixB: Baseline, optimal policy
    !v_matrixA: Alternative, suboptimal policy
    
    dev_v2 = v_matrixB(i_b0B, i_a0)-v_matrixA(i_b0A, i_a0)
    PRINT*,'Indices of conditioning state (a,b) =',i_a0,i_b0A,i_b0B
    PRINT*,'Cond values for Baseline, Alternative=',  v_matrixB(i_b0B, i_a0), v_matrixA(i_b0A, i_a0)
    PRINT*,'Initial V_baseline-V_alt =', dev_v2,gain_c
    dev_v2_old = dev_v2
    PRINT*,' '

do WHILE(convergence2<0)
    iter=iter+1
    iter2 = 0
    convergence=-1

    v0_matrixA_w  = v0_matrixA
    v_matrixA_w   = v_matrixA
    vaut_matrixA_w= vaut_matrixA
    
    do WHILE(convergence<0)     !convergence given cons gain lambda
        iter2 = iter2+1    

    do i_a = 1,na         
         IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
             iii_b=index_bfeas4(i_a)+1
             nbaux =nb-iii_b+1
             FDATA(1:nbaux) = v0_matrixA_w(iii_b:nb,i_a)
	         ILEFT  = 0
	         IRIGHT = 0
	         DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(iii_b+1) - bgrid(iii_b))
	         DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	         call  DCSDEC (nbaux, bgrid(iii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix_w(i_a, iii_b:nb), coeff_matrix_w(i_a, :,iii_b:nb))
         ENDIF         
     ENDDO

    dev_v = 0d+0
                
    do i_a = 1,na
        do i_b = 1,nb

            a_initial = agrid(i_a)
            b_initial = bgrid(i_b)
            
            nd2=2
            IF (psi1.GE.500d0) nd2=1
            
            do i_d =1,nd2
                i_default_today = i_d
    
                IF (i_default_today.EQ.1) THEN
                    cres  = (1d0+gain_c)*c_matrixA(i_b, i_a)
                    gnres = (1d0+gain_c)*gn_matrixA(i_b, i_a)
                    hres  = x_matrixA(i_b, i_a)
                    bpres = b_next_matrixA(i_b, i_a)
                ELSEIF (i_default_today.EQ.2) THEN
                    cres  = (1d0+gain_c)*caut_matrixA(i_b, i_a)
                    gnres = (1d0+gain_c)*gnaut_matrixA(i_b, i_a)
                    hres  = xaut_matrixA(i_b, i_a)
                    bpres = 0d0
                ENDIF
 
                aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
                aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma)) 
                
                IF (sigma.EQ.1d0) THEN
                    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
                ELSE    
                    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
                ENDIF                   
                    
                IF (sigmag.EQ.1d0) THEN 
                    utilg = DLOG(MAX(gnres,1d-8))
                ELSE
                    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
                ENDIF
    
                u_value = (1d0-chi)*utilc+chi*utilg

                a_threshold   = a_fun(bpres,a_initial)            !using original value fcns
                cdf_threshold = DNORDF((a_threshold - rhoa*a_initial - (1-rhoa)* mua)/sigmaa)

                exp_v_next = 0d0
                exp_v_next_aut = 0d0
                scalar = 1d0 / SQRT(pi_number)

                a_min1 = (1-rhoa)*mua + rhoa* a_initial - 4*sigmaa
                a_max1 = (1-rhoa)*mua + rhoa* a_initial + 4*sigmaa

                NINTV = ncdf - 1
                acum1=0d+0
                acum2=0d+0

               if (a_threshold < a_min1) then

                   do i=1,neps 
                      cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                      a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                      value_next = v0_fun_welfare(bpres,a_next)
                      acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
                   end do

               elseif(a_threshold > a_max1) then
                   do i=1,neps
                      cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                      a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                      value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                      acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                   end do

               else
                  do i=1,neps

                    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
                    a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                    value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next
                    
            !        compute integral for a > a_threshold
                    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
                    a_next1 = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf1) * sigmaa 
                    value_next = v0_fun_welfare(bpres,a_next1)
                    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next
                                
                  end do
               end if

            exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

            IF (i_default_today.EQ.2) THEN
                acum3=0d+0
                    do i=1,neps
                        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                        a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                        value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                        acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                    end do
    
                    exp_v_next_aut = acum3/(cdfmax - cdfmin)
                    exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
            ENDIF

            ucost = defaultgrid(i_default_today)*udef(a_initial)
            objective_fun = (u_value - ucost) + beta * exp_v_next

            IF (i_default_today.EQ.2) THEN
                vaut_matrixA_w_new(i_b,i_a) = objective_fun
            ELSE
                v0_matrixA_w_new(i_b,i_a) = objective_fun         
            ENDIF

            ENDDO       !i_d

            i_defaultA = default_decisionA(i_b, i_a)
            defval = defaultgrid(i_defaultA)
            v_matrixA_w_new(i_b, i_a) = defval*vaut_matrixA_w_new(i_b,i_a) + (1d0-defval)*v0_matrixA_w_new(i_b,i_a) 

            v_matrixA_w_new(i_b, i_a) = MAX(v_matrixA_w_new(i_b, i_a),-1000d0)
            dev_v =  max(ABS(v_matrixA_w_new(i_b, i_a) - v_matrixA_w(i_b, i_a)), dev_v) 

ENDDO       !i_a
    ENDDO       !i_b
           
    vaut_matrixA_w = vaut_matrixA_w_new
    v0_matrixA_w   = v0_matrixA_w_new
    v_matrixA_w    = v_matrixA_w_new                  
                    
    if (DABS(dev_gain).GT.0.1d0) THEN
        if (dev_v < 1d-2) then
            convergence =1
        end if
    else
        if (dev_v < criteria) then
            convergence =1
        end if
    endif

    ENDDO

    dev_v2_old = dev_v2  
    dev_v2 =  v_matrixB(i_b0B, i_a0) - v_matrixA_w(i_b0A, i_a0)

    if (DABS(dev_v2) < criteria) then
        convergence2 =1
        PRINT*,'ii_b, consumption gain',ii_b, gain_c
            gain_c_state(i_b0) = gain_c
                WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10) ')  bgrid(i_b0),  gain_c_state(i_b0),  v_matrixB(i_b0B, i_a0), v_matrixA(i_b0A, i_a0)

            PRINT*,'_________________________________________________'
            PRINT*,'b, Cond consumption gain =', ii_b,gain_c
            PRINT*,'_________________________________________________'
    end if

    i_defaultA = default_decisionA(i_b0, i_a0)
    defval = defaultgrid(i_defaultA)
    
    IF (defval.EQ.1d0) THEN
        c0     = (1d0+gain_C)*caut_matrixA(i_b0, i_a0)
        gn0    = (1d0+gain_C)*gnaut_matrixA(i_b0, i_a0)
        hres   = xaut_matrixA(i_b0, i_a0)
        aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
        aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))    
        ELSE
        c0     = (1d0+gain_C)*c_matrixA(i_b0, i_a0)       
        gn0    = (1d0+gain_C)*gn_matrixA(i_b0, i_a0)
        hres   = x_matrixA(i_b0, i_a0)
        aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
        aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))    
    ENDIF

    aa0    = (aahat**(1d0-sigma))*aa1
    slope = (1d0-chi)*aa0*c0**(1d0-sigma)+chi*gn0**(1d0-sigma)
    slope = slope/(1d0-beta)

    gain_c_new = (v_matrixB(i_b0B, i_a0) - v_matrixA_w(i_b0A, i_a0))/slope + gain_c

    dev_gain = gain_c_new -gain_c
    weight_gain=DABS(dev_v2_old)/(DABS(dev_v2_old)+DABS(dev_v2))

    WRITE(*, '(I5, X, F12.8, X, F12.8, X, F12.8, X, F12.8, X, F12.8)') iter, gain_c, v_matrixA_w(i_b0A, i_a0), dev_gain,dev_v2,slope

    lag_lim=-0.9d0
    gain_c   = weight_gain*MAX(gain_c_new,lag_lim)+(1d0-weight_gain)*gain_c

ENDDO
ENDDO

CLOSE(1111)

PRINT*,'WELFARE GAINS COMPUTED'


END SUBROUTINE WELFARE_GAIN
!******************************************************************
SUBROUTINE WELFARE_GAIN_SAM
!Computes welfare gains in terms of permanent increase in total consumption from switching from 
!Samuelson rule to Baseline, as function of current yT
! Economy A: to which add welfare gain (Samuelson)
! Economy B: Better economy (Baseline)

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

REAL(8),DIMENSION(nb,na)::v_matrixA,v0_matrixA,vaut_matrixA
REAL(8),DIMENSION(nb,na)::b_next_matrixA,q_matrixA,q_matrix_nodefA
REAL(8),DIMENSION(nb,na)::x_matrixA,gn_matrixA,xaut_matrixA,gnaut_matrixA
INTEGER,DIMENSION(nb,na)::default_decisionA

REAL(8),DIMENSION(nb,na)::cn_matrixA,cnaut_matrixA,ct_matrixA,ctaut_matrixA,c_matrixA,caut_matrixA
REAL(8),DIMENSION(nb,na)::tau_matrixA,tauaut_matrixA

REAL(8),DIMENSION(61,na)::v_matrixB,v0_matrixB,vaut_matrixB
REAL(8),DIMENSION(nb,na)::vaut_matrixA_w,vaut_matrixA_w_new,v0_matrixA_w_new,v_matrixA_w,v_matrixA_w_new
REAL(8)::dev_v,dev_v2,ctval,cnval,b_next,bpres,p0,q,mstar,a0,gain_c,criteria
REAL(8)::hres,cres,gnres,taures,utilc,utilg,u_value,a_initial,b_initial,aux0

REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::exp_vA,exp_vB,a00,b00 
INTEGER,PARAMETER::nbA=nb2,nbB=nb2
REAL(8)::distAaccess(nbA,na),distBaccess(nbB,na)
REAL(8),DIMENSION(na)::distAaut,distBaut

REAL(8)::a_threshold,a_max1,a_min1,a_next,a_next1,acum2,acum1,acum3,cdf,cdf1,cdf_threshold,scalar
REAL(8)::exp_v_next,exp_v_next_aut,a_fun,value_next,ucost,udef,DNORDF,objective_fun,defval
REAL(8)::c0,gn0,slope,DNORIN
EXTERNAL::a_fun,udef

INTEGER ILEFT, IRIGHT,NINTV,choice_cond,i_a_global,i_b_global,iter,iter2,nbtemp
EXTERNAL DNORIN
INTEGER::i,i_a,i_b,ii_b,iii_b,i_d,i0(1),i_a0,i_b0,i_b0A,i_b0B,i_default_today,i_defaultA,convergence,convergence2
REAL(8)::DCSVAL,gain_c_new,dev_gain ,valuex,dev_v2_old
EXTERNAL::DCSVAL

REAL(8)::v0_fun,gain_c_state(nb),v0_fun_welfare,weight_gain,bgridB(61)
INTEGER::nbaux,index_opt,nd2

EXTERNAL::v0_fun,v0_fun_welfare
REAL(8)::new_value,hval,c01,c02,aahat,aa1,aa0,a01,a02 ,slope1,slope2,lag_lim 
REAL(8)::exp_vA_repay,exp_vB_repay ,v22, valueA,valueB, slope_aux(nb)


open (110, FILE='graphs\b_grid.txt')
DO i=1,nb
    READ(110,*) bgrid(i)
ENDDO

CLOSE(110)

!NOTE: To compute policy fcns with (na,nb), need to reset nb2, re-run func_gN first and store txt files for policy fcns in folder welfare
open (210, FILE='graphs_baselineX\v.txt')
open (211, FILE='graphs_baselineX\bgrid.txt')
do i_b = 1,61
    READ(211,*) bgridB(i_b)
    do i_a = 1,na
        READ(210, '(F18.10, X, F18.10, X, F18.10)') v_matrixB(i_b, i_a), &
                    v0_matrixB(i_b, i_a), vaut_matrixB(i_b, i_a)
    end do
end do

CLOSE(210)
CLOSE(211)

open (10, FILE='graphs\v.txt')
open (11, FILE='graphs\default.txt')
open (13, FILE='graphs\b_next.txt')
open (114, FILE='output\welfare\gh.txt')
open (115, FILE='output\welfare\ggn.txt')
open (116, FILE='output\welfare\ghaut.txt')
open (117, FILE='output\welfare\ggnaut.txt')
open (118, FILE='output\welfare\gct.txt')
open (119, FILE='output\welfare\gctaut.txt')

do i_a = 1,na
do i_b = 1,nb
    READ(114, *) x_matrixA(i_b, i_a)
    READ(115, *) gn_matrixA(i_b, i_a)                     
    READ(118, *) ct_matrixA(i_b, i_a)

ENDDO
ENDDO



do i_b = 1,nb
do i_a = 1,na
    i_a_global = i_a
    i_b_global = i_b

    READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrixA(i_b, i_a), &
                v0_matrixA(i_b, i_a), vaut_matrixA(i_b, i_a)
    READ(13, '(F15.11)') b_next_matrixA(i_b, i_a)
    if (vaut_matrixA(i_b, i_a) > v0_matrixA(i_b, i_a)) then
        default_decisionA(i_b, i_a) =2
    else
        default_decisionA(i_b, i_a) =1
    end if
end do
end do

do i_a = 1,na
    READ(116, *) valuex
    xaut_matrixA(:, i_a)=valuex
    READ(117, *) valuex
    gnaut_matrixA(:, i_a)=valuex 
    READ(119, *) valuex
    ctaut_matrixA(:, i_a)=valuex 
enddo    

CLOSE(10)
CLOSE(13)
CLOSE(114)
CLOSE(115)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)

PRINT*,' '
PRINT*,'SOLUTION UPLOADED FOR WELFARE ANALYSIS'
PRINT*,' '

!upload value and policy fcns
PRINT*,'Compute CONDITIONAL consumption gain' 
PRINT*,' '

do i_a = 1,na
do i_b = 1,nb
    a0 = DEXP(agrid(i_a))
    a_initial = agrid(i_a)
    b_initial = bgrid(i_b)
    
    hval = x_matrixA(i_b, i_a)
    cnval = hval**alpha-gn_matrixA(i_b, i_a)
    cn_matrixA(i_b, i_a)= cnval
    
    b_next = b_next_matrixA(i_b, i_a)
    ctval = ct_matrixA(i_b, i_a) 

    IF (mu .NE. 0d0) THEN
        c_matrixA(i_b, i_a) = (omegac*ct_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cn_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        c_matrixA(i_b, i_a) = (ct_matrixA(i_b, i_a)**(omegac))*(cn_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
    
    cnval = xaut_matrixA(i_b, i_a)**alpha-gnaut_matrixA(i_b, i_a)
    cnaut_matrixA(i_b, i_a)= cnval
    ctval = ctaut_matrixA(i_b, i_a)

    IF (mu .NE. 0d0) THEN
        caut_matrixA(i_b, i_a) = (omegac*ctaut_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cnaut_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        caut_matrixA(i_b, i_a) = (ctaut_matrixA(i_b, i_a)**(omegac))*(cnaut_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
ENDDO
ENDDO

OPEN(1111,FILE='output//gain_samuelson_yT.txt')

!To apply a_fun, need original value fcns
!Default thresholds computed using model A
vaut_matrix  = vaut_matrixA
v0_matrix    = v0_matrixA
v_matrix     = v_matrixA

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
    ii_b=index_bfeas4(i_a)+1
    nbaux =nb-ii_b+1
    FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
ENDIF         
ENDDO

criteria=1d-5
iter=0
gain_c_state = 0.01d0

a00 = mua/(1d0-rhoa)       !start from unconditional mean
i_a0 = LOCATE2(agrid,a00)

DO i_a0=1,na

i0 = LOCATE2(bgrid,0d0)
i_b0A = i0(1)
i0 = LOCATE2(bgridB,0d0)
i_b0B = i0(1)

DO ii_b=i_b0A,i_b0A
    i_b0 = ii_b

    IF (MIN(v_matrixB(i_b0B, i_a0),v_matrixA(i_b0A, i_a0)).LE.-500d0) THEN
        gain_c_state(i_b0) = -1d0
        CYCLE
    ENDIF

    convergence = -1
    convergence2= -1

    gain_c = 0.01d0
    dev_gain  = 1d0

    !v_matrixB: Baseline, optimal policy
    !v_matrixA: Alternative, suboptimal policy
    
    dev_v2 = v_matrixB(i_b0B, i_a0)-v_matrixA(i_b0A, i_a0)
    PRINT*,'Indices of conditioning state (a,b) =',i_a0,i_b0,i_b0B
    PRINT*,'Cond values for Baseline, Alternative=',  v_matrixB(i_b0B, i_a0), v_matrixA(i_b0A, i_a0)
    PRINT*,'Initial V_baseline-V_alt =', dev_v2,gain_c
    dev_v2_old = dev_v2
    PRINT*,' '

do WHILE(convergence2<0)
    iter=iter+1
    iter2 = 0
    convergence=-1

    v0_matrixA_w  = v0_matrixA
    v_matrixA_w   = v_matrixA
    vaut_matrixA_w= vaut_matrixA
    
    do WHILE(convergence<0)     !convergence given cons gain lambda
        iter2 = iter2+1    

    do i_a = 1,na         
         IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
             iii_b=index_bfeas4(i_a)+1
             nbaux =nb-iii_b+1
             FDATA(1:nbaux) = v0_matrixA_w(iii_b:nb,i_a)
	         ILEFT  = 0
	         IRIGHT = 0
	         DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(iii_b+1) - bgrid(iii_b))
	         DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	         call  DCSDEC (nbaux, bgrid(iii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix_w(i_a, iii_b:nb), coeff_matrix_w(i_a, :,iii_b:nb))
         ENDIF         
     ENDDO

    dev_v = 0d+0
                
    do i_a = 1,na
        do i_b = 1,nb

            a_initial = agrid(i_a)
            b_initial = bgrid(i_b)
            
            nd2=2
            IF (psi1.GE.500d0) nd2=1
            
            do i_d =1,nd2
                i_default_today = i_d
    
                IF (i_default_today.EQ.1) THEN
                    cres  = (1d0+gain_c)*c_matrixA(i_b, i_a)
                    gnres = (1d0+gain_c)*gn_matrixA(i_b, i_a)
                    hres  = x_matrixA(i_b, i_a)
                    bpres = b_next_matrixA(i_b, i_a)
                ELSEIF (i_default_today.EQ.2) THEN
                    cres  = (1d0+gain_c)*caut_matrixA(i_b, i_a)
                    gnres = (1d0+gain_c)*gnaut_matrixA(i_b, i_a)
                    hres  = xaut_matrixA(i_b, i_a)
                    bpres = 0d0
                ENDIF
 
                aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
                aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma)) 
                
                IF (sigma.EQ.1d0) THEN
                    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
                ELSE    
                    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
                ENDIF                   
                    
                IF (sigmag.EQ.1d0) THEN 
                    utilg = DLOG(MAX(gnres,1d-8))
                ELSE
                    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
                ENDIF
    
                u_value = (1d0-chi)*utilc+chi*utilg

                a_threshold   = a_fun(bpres,a_initial)            !using original value fcns
                cdf_threshold = DNORDF((a_threshold - rhoa*a_initial - (1-rhoa)* mua)/sigmaa)

                exp_v_next = 0d0
                exp_v_next_aut = 0d0
                scalar = 1d0 / SQRT(pi_number)

                a_min1 = (1-rhoa)*mua + rhoa* a_initial - 4*sigmaa
                a_max1 = (1-rhoa)*mua + rhoa* a_initial + 4*sigmaa

                NINTV = ncdf - 1
                acum1=0d+0
                acum2=0d+0

               if (a_threshold < a_min1) then

                   do i=1,neps 
                      cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                      a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                      value_next = v0_fun_welfare(bpres,a_next)
                      acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
                   end do

               elseif(a_threshold > a_max1) then
                   do i=1,neps
                      cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                      a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                      value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                      acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                   end do

               else
                  do i=1,neps

                    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
                    a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                    value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next
                    
            !        compute integral for a > a_threshold
                    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
                    a_next1 = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf1) * sigmaa 
                    value_next = v0_fun_welfare(bpres,a_next1)
                    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next
                                
                  end do
               end if

            exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

            IF (i_default_today.EQ.2) THEN
                acum3=0d+0
                    do i=1,neps
                        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                        a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                        value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                        acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                    end do
    
                    exp_v_next_aut = acum3/(cdfmax - cdfmin)
                    exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
            ENDIF

            ucost = defaultgrid(i_default_today)*udef(a_initial)
            objective_fun = (u_value - ucost) + beta * exp_v_next

            IF (i_default_today.EQ.2) THEN
                vaut_matrixA_w_new(i_b,i_a) = objective_fun
            ELSE
                v0_matrixA_w_new(i_b,i_a) = objective_fun         
            ENDIF

            ENDDO       !i_d

            i_defaultA = default_decisionA(i_b, i_a)
            defval = defaultgrid(i_defaultA)
            v_matrixA_w_new(i_b, i_a) = defval*vaut_matrixA_w_new(i_b,i_a) + (1d0-defval)*v0_matrixA_w_new(i_b,i_a) 

            v_matrixA_w_new(i_b, i_a) = MAX(v_matrixA_w_new(i_b, i_a),-1000d0)
            dev_v =  max(ABS(v_matrixA_w_new(i_b, i_a) - v_matrixA_w(i_b, i_a)), dev_v) 

ENDDO       !i_a
    ENDDO       !i_b
           
    vaut_matrixA_w = vaut_matrixA_w_new
    v0_matrixA_w   = v0_matrixA_w_new
    v_matrixA_w    = v_matrixA_w_new                  
                    
    if (DABS(dev_gain).GT.0.1d0) THEN
        if (dev_v < 1d-2) then
            convergence =1
        end if
    else
        if (dev_v < criteria) then
            convergence =1
        end if
    endif

    ENDDO

    dev_v2_old = dev_v2  
    dev_v2 =  v_matrixB(i_b0B, i_a0) - v_matrixA_w(i_b0A, i_a0)

    if (DABS(dev_v2) < criteria) then
        convergence2 =1
        PRINT*,'ii_b, consumption gain',ii_b, gain_c
            gain_c_state(i_b0) = gain_c
                WRITE(1111,' (F18.10, X, F18.10, X, F18.10, X, F18.10) ')  bgrid(i_b0),  gain_c_state(i_b0),  v_matrixB(i_b0B, i_a0), v_matrixA(i_b0A, i_a0)

            PRINT*,'_________________________________________________'
            PRINT*,'b, Cond consumption gain =', ii_b,gain_c
            PRINT*,'_________________________________________________'
    end if

    i_defaultA = default_decisionA(i_b0, i_a0)
    defval = defaultgrid(i_defaultA)
    
    IF (defval.EQ.1d0) THEN
        c0     = (1d0+gain_C)*caut_matrixA(i_b0, i_a0)
        gn0    = (1d0+gain_C)*gnaut_matrixA(i_b0, i_a0)
        hres   = xaut_matrixA(i_b0, i_a0)
        aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
        aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))    
        ELSE
        c0     = (1d0+gain_C)*c_matrixA(i_b0, i_a0)       
        gn0    = (1d0+gain_C)*gn_matrixA(i_b0, i_a0)
        hres   = x_matrixA(i_b0, i_a0)
        aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
        aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))    
    ENDIF

    aa0    = (aahat**(1d0-sigma))*aa1
    slope = (1d0-chi)*aa0*c0**(1d0-sigma)+chi*gn0**(1d0-sigma)
    slope = slope/(1d0-beta)

    IF (slope.LT.0d0) THEN
        PRINT*,'Problem slope ',slope
        PRINT*,'slope1,slope2',slope1,slope2           
        PRINT*,'a01,c01',aa0,c0
        PRINT*,'defval,',defaultgrid(default_decisionA(i_b0, i_a0)),defaultgrid(default_decisionA(1,1)),defaultgrid(default_decisionA(nb,na))
        PRINT*,'cc',c_matrixA(i_b0, i_a0),caut_matrixA(i_b0, i_a0)
        PRINT*,i_b0, i_a0
        PAUSE
        STOP
    ENDIF        
    
    gain_c_new = (v_matrixB(i_b0B, i_a0) - v_matrixA_w(i_b0A, i_a0))/slope + gain_c

    dev_gain = gain_c_new -gain_c
    weight_gain=DABS(dev_v2_old)/(DABS(dev_v2_old)+DABS(dev_v2))

    WRITE(*, '(I5, X, I5, X, F12.8, X, F12.8, X, F12.8, X, F12.8)') ii_b,iter, gain_c, dev_gain,dev_v2,slope
    WRITE(1111, '(I5,  X, F12.8, X, F12.8)') i_a0,agrid(i_a0), gain_c

    lag_lim=-0.9d0
    gain_c   = weight_gain*MAX(gain_c_new,lag_lim)+(1d0-weight_gain)*gain_c

ENDDO
ENDDO
ENDDO

CLOSE(1111)

PRINT*,'WELFARE GAINS COMPUTED'


END SUBROUTINE WELFARE_GAIN_SAM
!******************************************************************
SUBROUTINE WELFARE_GAIN_UNCOND
! Economy A: to which add welfare gain (Samuelson)
! Economy B: Better economy (Baseline)


USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

INTEGER,PARAMETER::nbA=nb2,nbB=nb2
REAL(8),DIMENSION(nb,na)::v_matrixA,v0_matrixA,vaut_matrixA
REAL(8),DIMENSION(nb,na)::x_matrixA,gn_matrixA,xaut_matrixA,gnaut_matrixA,b_next_matrixA
REAL(8),DIMENSION(nb,na)::cn_matrixA,cnaut_matrixA,ct_matrixA,ctaut_matrixA,c_matrixA,caut_matrixA
INTEGER,DIMENSION(nb,na)::default_decisionA

REAL(8),DIMENSION(nbB,na)::v_matrixB,v0_matrixB(nbB, na),vaut_matrixB(nbB,na) 
REAL(8),DIMENSION(nb,na)::vaut_matrixA_w,vaut_matrixA_w_new,v0_matrixA_w_new,v_matrixA_w,v_matrixA_w_new
REAL(8)::dev_v,dev_v2,ctval,cnval,b_next,bpres,a0,gain_c,criteria
REAL(8)::hres,cres,gnres,utilc,utilg,u_value,a_initial,b_initial

REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb)
REAL(8)::exp_vA,exp_vB,a00,b00 
REAL(8)::distAaccess(nbA,na),distBaccess(nbB,na)
REAL(8),DIMENSION(na)::distAaut,distBaut

REAL(8)::a_threshold,a_max1,a_min1,a_next,a_next1,acum2,acum1,acum3,cdf,cdf1,cdf_threshold,scalar
REAL(8)::exp_v_next,exp_v_next_aut,a_fun,value_next,ucost,udef,DNORDF,objective_fun,defval
REAL(8)::c0,gn0,slope,DNORIN,value1,value2
EXTERNAL::a_fun,udef

INTEGER ILEFT, IRIGHT,NINTV,iter,iter2,nbtemp
EXTERNAL DNORIN
INTEGER::i,i_a,i_b,ii_b,iii_b,i_d,i_a0,i_b0,i_default_today,i_defaultA,convergence,convergence2,nd2
REAL(8)::DCSVAL,gain_c_new,dev_gain ,valuex,dev_v2_old
EXTERNAL::DCSVAL

REAL(8)::v0_fun,v0_fun_welfare,weight_gain
INTEGER::nbaux,index_opt

EXTERNAL::v0_fun,v0_fun_welfare
REAL(8)::new_value,hval,c01,c02,aahat,aa1,aa0,a01,a02 ,slope1,slope2,lag_lim 
REAL(8)::exp_vA_repay,exp_vB_repay ,v22, valueA,valueB, slope_aux(nb)

open (210, FILE='output_baselineX\welfare_uncond\gvcon.txt')
open (211, FILE='output_baselineX\welfare_uncond\gvaut.txt')

do i_a = 1,na
    READ(211,*) value1
    vaut_matrixB(:, i_a) = value1
    do i_b = 1,nbB
    READ(210,*) value2
    v0_matrixB(i_b, i_a)=value2
    v_matrixB(i_b, i_a)=MAX(value1,value2)
ENDDO
ENDDO

CLOSE(210)
CLOSE(211)

OPEN(1,FILE='output\bgrid2.txt')    
            
DO i=1,nb2
    READ(1,*) bgrid2(i)
ENDDO

CLOSE(1)

open (10, FILE='graphs\v.txt')
open (11, FILE='graphs\default.txt')
open (13, FILE='graphs\b_next.txt')
open (114, FILE='output\welfare\gh.txt')
open (115, FILE='output\welfare\ggn.txt')
open (116, FILE='output\welfare\ghaut.txt')
open (117, FILE='output\welfare\ggnaut.txt')
open (118, FILE='output\welfare\gct.txt')
open (119, FILE='output\welfare\gctaut.txt')

do i_a = 1,na
do i_b = 1,nb
    READ(114, *) x_matrixA(i_b, i_a)
    READ(115, *) gn_matrixA(i_b, i_a)
    READ(118, *) ct_matrixA(i_b, i_a)
ENDDO

ENDDO
    
do i_b = 1,nb
do i_a = 1,na
    READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrixA(i_b, i_a), &
                v0_matrixA(i_b, i_a), vaut_matrixA(i_b, i_a)
    READ(13, '(F15.11)') b_next_matrixA(i_b, i_a)
    if (vaut_matrixA(i_b, i_a) > v0_matrixA(i_b, i_a)) then
        default_decisionA(i_b, i_a) =2
    else
        default_decisionA(i_b, i_a) =1
    end if
end do
end do

do i_a = 1,na
    READ(116, '(F15.11)') valuex
    xaut_matrixA(:, i_a)=valuex
    READ(117, '(F15.11)') valuex
    gnaut_matrixA(:, i_a)=valuex 
    READ(119, '(F15.11)') valuex
    ctaut_matrixA(:, i_a)=valuex 
enddo    

CLOSE(10)
CLOSE(13)
CLOSE(114)
CLOSE(115)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)

PRINT*,' '
PRINT*,'SOLUTION UPLOADED FOR WELFARE ANALYSIS'
PRINT*,' '

!upload value and policy fcns
nbtemp = nb

CALL HIST(distAaccess,distBaccess,distAaut,distBaut)
distBaccess=distAaccess !use flex-wage gn distribution to compute expected values

nbtemp = 1
PRINT*,'Compute UNCONDITIONAL consumption gain' 

do i_a = 1,na
do i_b = 1,nb
    a0 = DEXP(agrid(i_a))
    a_initial = agrid(i_a)
    b_initial = bgrid(i_b)
    
    hval = x_matrixA(i_b, i_a)
    cnval = hval**alpha-gn_matrixA(i_b, i_a)
    cn_matrixA(i_b, i_a)= cnval
    b_next = b_next_matrixA(i_b, i_a)
    ctval = ct_matrixA(i_b, i_a) 

    IF (mu .NE. 0d0) THEN
        c_matrixA(i_b, i_a) = (omegac*ct_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cn_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        c_matrixA(i_b, i_a) = (ct_matrixA(i_b, i_a)**(omegac))*(cn_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
    
    cnval = xaut_matrixA(i_b, i_a)**alpha-gnaut_matrixA(i_b, i_a)
    cnaut_matrixA(i_b, i_a)= cnval
    ctval = ctaut_matrixA(i_b, i_a)

    IF (mu .NE. 0d0) THEN
        caut_matrixA(i_b, i_a) = (omegac*ctaut_matrixA(i_b, i_a)**(-mu)+ (1d0-omegac)*cnaut_matrixA(i_b, i_a)**(-mu))**(-1d0/mu)
    ELSE
        caut_matrixA(i_b, i_a) = (ctaut_matrixA(i_b, i_a)**(omegac))*(cnaut_matrixA(i_b, i_a)**(1d0-omegac))
    ENDIF
ENDDO
ENDDO

!To apply a_fun, need original value fcns
!Default thresholds computed using model A
vaut_matrix  = vaut_matrixA
v0_matrix    = v0_matrixA
v_matrix     = v_matrixA

index_bfeas4 = 0d0
    
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
    ii_b=index_bfeas4(i_a)+1
    nbaux =nb-ii_b+1
    FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
ENDIF         
ENDDO

criteria=1d-5
iter=0

a00 = mua/(1d0-rhoa) !start from unconditional mean
i_a0 = LOCATE2(agrid,a00)


b00  = -1.d0*1.108d0*(lambda+r) !condition on asymptotic debt
i_b0 = LOCATE2(bgrid,b00)

convergence = -1
convergence2= -1

gain_c = 0.02826d0
dev_gain  = 1d0
    
exp_vA = 0d0
exp_vB = 0d0
slope  = 0d0
    
!to compute expected values use debt grids from policy fcns
DO i_a=1,na2
    DO i_b=1,nb2
        valueA=LINEAR_INT(bgrid,v0_matrixA(:, i_a),bgrid2(i_b))
!        valueB=LINEAR_INT(bgrid,v0_matrixB(:, i_a),bgrid2(i_b))
        valueB=v0_matrixB(i_b, i_a)

        IF (distAaccess(i_b,i_a).GT.1d-8) exp_vA = exp_VA + distAaccess(i_b,i_a)*valueA
        IF (distBaccess(i_b,i_a).GT.1d-8) exp_vB = exp_VB + distBaccess(i_b,i_a)*valueB
        PRINT*,i_b,i_a,valueA,valueB
        PRINT*,exp_vA ,exp_vB ,distAaccess(i_b,i_a),distBaccess(i_b,i_a)
    ENDDO

    exp_vA = exp_VA + distAaut(i_a)*vaut_matrixA(1, i_a)
    exp_vB = exp_VB + distBaut(i_a)*vaut_matrixB(1, i_a)
    PRINT*,'A',i_a,vaut_matrixA(1, i_a),vaut_matrixB(1, i_a)
    PRINT*,exp_vA ,exp_vB ,distAaut(i_a),distBaut(i_a)
        
ENDDO
    
dev_v2 =  exp_VB-exp_VA
PRINT*,'Unconditional Baseline, Alternative=', exp_VB, exp_VA
PRINT*,'Initial difference =', dev_v2

dev_v2_old = dev_v2
PRINT*,' '
    
do WHILE(convergence2<0)
    iter=iter+1
    iter2 = 0
    convergence=-1

    v0_matrixA_w  = v0_matrixA
    v_matrixA_w   = v_matrixA
    vaut_matrixA_w= vaut_matrixA            

    PRINT*,'gain_c',gain_c
    
    do WHILE(convergence<0)     !convergence given cons gain lambda
        iter2 = iter2+1    

        do i_a = 1,na         
             IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
                 iii_b=index_bfeas4(i_a)+1
                 nbaux =nb-iii_b+1
                 FDATA(1:nbaux) = v0_matrixA_w(iii_b:nb,i_a)
	             ILEFT  = 0
	             IRIGHT = 0
	             DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(iii_b+1) - bgrid(iii_b))
	             DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	             call  DCSDEC (nbaux, bgrid(iii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix_w(i_a, iii_b:nb), coeff_matrix_w(i_a, :,iii_b:nb))
             ENDIF         
         ENDDO

        dev_v = 0d+0
                
        do i_a = 1,na
            do i_b = 1,nb

                a_initial = agrid(i_a)
                b_initial = bgrid(i_b)
            
                nd2=2
                IF (psi1.GE.500d0) nd2=1

                do i_d =1,nd2
                    i_default_today = i_d
    
                    IF (i_default_today.EQ.1) THEN
                        cres  = (1d0+gain_c)*c_matrixA(i_b, i_a)
                        gnres = (1d0+gain_c)*gn_matrixA(i_b, i_a)
                        hres  = x_matrixA(i_b, i_a)
                        bpres = b_next_matrixA(i_b, i_a)
                    ELSEIF (i_default_today.EQ.2) THEN
                        cres  = (1d0+gain_c)*caut_matrixA(i_b, i_a)
                        gnres = (1d0+gain_c)*gnaut_matrixA(i_b, i_a)
                        hres  = xaut_matrixA(i_b, i_a)
                        bpres = 0d0
                    ENDIF
 
                    aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
                    aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))
                
                    IF (sigma.EQ.1d0) THEN
                        utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
                    ELSE    
                        utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
                    ENDIF                   
                    
                    IF (sigmag.EQ.1d0) THEN 
                        utilg = DLOG(MAX(gnres,1d-8))
                    ELSE
                        utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
                    ENDIF
    
                    u_value = (1d0-chi)*utilc+chi*utilg
    
                    IF (psi1.LT.500d0) THEN
                        a_threshold   = a_fun(bpres,a_initial)            !using original value fcns
                    ELSE
                        a_threshold   = -1d4
                    ENDIF

                    cdf_threshold = DNORDF((a_threshold - rhoa*a_initial - (1-rhoa)* mua)/sigmaa)

                    exp_v_next = 0d0
                    exp_v_next_aut = 0d0
                    scalar = 1d0 / SQRT(pi_number)

                    a_min1 = (1-rhoa)*mua + rhoa* a_initial - 4*sigmaa
                    a_max1 = (1-rhoa)*mua + rhoa* a_initial + 4*sigmaa

                    NINTV = ncdf - 1
                    acum1=0d+0
                    acum2=0d+0

                   if (a_threshold < a_min1) then

                       do i=1,neps 
                          cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                          a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                          value_next = v0_fun_welfare(bpres,a_next)
                          acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
                       end do

                   elseif(a_threshold > a_max1) then
                       do i=1,neps
                          cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                          a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                          value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                          acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                       end do

                   else
                       !v22=0d0
                      do i=1,neps

                        cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
                        a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                        value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                        acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next
                        !v22=v22+0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next
                    
                !        compute integral for a > a_threshold
                        cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
                        a_next1 = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf1) * sigmaa 
                        value_next = v0_fun_welfare(bpres,a_next1)
                        acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next
                                
                      end do
                   end if

                exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

                IF (i_default_today.EQ.2) THEN
                    acum3=0d+0
                        do i=1,neps
                            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
                            a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
                            value_next = LINEAR_INT(agrid,vaut_matrixA_w(1,:),a_next)
                            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
                        end do
    
                        exp_v_next_aut = acum3/(cdfmax - cdfmin)
                        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
                ENDIF

                ucost = defaultgrid(i_default_today)*udef(a_initial)

                objective_fun = (u_value - ucost) + beta * exp_v_next
                                    
                IF (i_default_today.EQ.2) THEN
                    vaut_matrixA_w_new(i_b,i_a) = objective_fun
                ELSE
                    v0_matrixA_w_new(i_b,i_a) = objective_fun
                ENDIF

            ENDDO

            i_defaultA = default_decisionA(i_b, i_a)
            defval = defaultgrid(i_defaultA)
            v_matrixA_w_new(i_b, i_a) = defval*vaut_matrixA_w_new(i_b,i_a) + (1d0-defval)*v0_matrixA_w_new(i_b,i_a) 

            dev_v =  max(ABS(v_matrixA_w_new(i_b, i_a) - v_matrixA_w(i_b, i_a)), dev_v)
        
        ENDDO       !i_a
    ENDDO       !i_b
           
    vaut_matrixA_w = vaut_matrixA_w_new
    v0_matrixA_w   = v0_matrixA_w_new
    v_matrixA_w    = v_matrixA_w_new     
    
    
    if (DABS(dev_gain).GT.0.1d0) THEN
        if (dev_v < 1d-2) then
            convergence =1
        end if
    else
        
        if (dev_v < criteria) then
        convergence =1
        end if
    endif

    WRITE(*, '(I5, X, F12.6, X, F12.6, X, F13.6)') iter2, gain_c,dev_v,dev_gain

    ENDDO
    
        exp_vA = 0d0
        slope  = 0d0        !proxy for slope: E u'(c)/(1-beta) - use representative-agent consumption

        DO i_a=1,na
        DO i_b=1,nb2
            valueA=LINEAR_INT(bgrid,v0_matrixA_w(:, i_a),bgrid2(i_b))
            exp_vA = exp_VA + distAaccess(i_b,i_a)*valueA
        ENDDO
        ENDDO
    
        slope_aux = 0d0
    
        DO i_a=1,na
        DO i_b=1,nb
            c0 = (1d0+gain_c)*c_matrixA(i_b, i_a)
            c0 = MAX(1d-12,c0)
            gn0 = (1d0+gain_c)*gn_matrixA(i_b, i_a)
            gn0 = MAX(1d-12,gn0)
            hres   = x_matrixA(i_b, i_a)
            aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
            aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))
        
            slope_aux(i_b) = ((1d0-chi)*(aahat**(1d0-sigma))*aa1*c0**(1d0-sigma)+chi*gn0**(1d0-sigma))
            
            IF (slope.LT.0d0) THEN
                PRINT*,'Problem slope Repay',slope
                PRINT*,'aahat,aa1,c0',(aahat**(1d0-sigma)),aa1,c0**(-sigma)*c0
                PAUSE
                STOP
            ENDIF
            ENDDO

            DO i_b=1,nb2
                slope = slope + distAaccess(i_b,i_a)*LINEAR_INT(bgrid,slope_aux,bgrid2(i_b))
            ENDDO

            exp_vA = exp_VA + distAaut(i_a)*vaut_matrixA_w(1, i_a)
            c0 = (1d0+gain_c)*caut_matrixA(1, i_a)
            gn0 = (1d0+gain_c)*gnaut_matrixA(1, i_a)
            hres   = xaut_matrixA(1, i_a)
            aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
            aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))

            slope = slope + distAaut(i_a)*((1d0-chi)*(aahat**(1d0-sigma))*aa1*c0**(1d0-sigma)+chi*gn0**(1d0-sigma))
        
            IF (slope.LT.0d0) THEN
                PRINT*,'Problem slope Aut',slope
                PRINT*,'aahat,aa1,c0',(aahat**(1d0-sigma)),aa1,c0**(-sigma)*c0
                PAUSE
                STOP
            ENDIF                   
        ENDDO

        slope = slope/(1d0-beta)
        dev_v2_old = dev_v2  
        dev_v2     =  exp_VB-exp_VA


if (DABS(dev_v2) < criteria) then
    convergence2 =1
    PRINT*,'ii_b, consumption gain',ii_b, gain_c
    PRINT*,'_________________________________________________'
    PRINT*,'Uncond consumption gain =', gain_c
    PRINT*,'_________________________________________________'
    RETURN
ELSE
    gain_c_new = (exp_VB-exp_VA)/slope + gain_c

    dev_gain = gain_c_new -gain_c
    weight_gain=DABS(dev_v2_old)/(DABS(dev_v2_old)+DABS(dev_v2))

    lag_lim=-0.9d0
    gain_c   = weight_gain*MAX(gain_c_new,lag_lim)+(1d0-weight_gain)*gain_c
    WRITE(*, '(F12.6, X, F12.6, X, F12.6, X, F12.6)') gain_c,(exp_VB-exp_VA),exp_VB,exp_VA
ENDIF


ENDDO

END SUBROUTINE WELFARE_GAIN_UNCOND
!*********************************************************************************
DOUBLE PRECISION function v0_fun_welfare(b,a)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: b,a
DOUBLE PRECISION ::  slope, DCSVAL, value_left, value_right 
DOUBLE PRECISION::bvec(3),avec(3),lamb(3),value1,value2,value3,interpolate

INTEGER :: index_a,index_b,ii_b,index_bfeas,index_bfeas3,NINTV,val,nbaux,index_bfeasmax,index1,index2
EXTERNAL DCSVAL,interpolate

!uncomment to use linear interpolation
!v0_fun=interpolate(b,a,v0_matrix)
!RETURN

NINTV = nb - 1
index_a = (MAX(MIN(INT((na-1)*(a-agrid(1))/(agrid(na)-agrid(1)))+1,na-1), 1))       !closest lower grid point

!realization of a belonds to agrid
IF (a-agrid(index_a).EQ.0d0) THEN
    index1=index_bfeas4(index_a)
    val=1
    IF (((b.GE.bgrid(index1+1)).AND.(index1+1.LT.nb)).OR.(index1.EQ.0)) THEN
    
        ii_b=index1+1
        nbaux =nb-ii_b+1
        v0_fun_welfare  = DCSVAL(b, nbaux-1, break_matrix_w(index_a,ii_b:nb), coeff_matrix_w(index_a,:,ii_b:nb))
        RETURN
    ELSEIF ((b.GE.bgrid(index1+1)).AND.(index1+1.EQ.nb)) THEN
        v0_fun_welfare  = v0_matrixA_w(nb,index_a)
    ELSE
        v0_fun_welfare =-1d6 
    ENDIF
    RETURN
ENDIF

!realization of a does not belong to agrid

index1 = index_bfeas4(index_a)
index2 = index_bfeas4(index_a+1)


IF (index1.EQ.0) THEN
    if (index2.EQ.0) THEN
        ii_b=index1+1
        nbaux =nb-ii_b+1
        
        value_left  = DCSVAL(b, nbaux-1, break_matrix_w(index_a,ii_b:nb), coeff_matrix_w(index_a,:,ii_b:nb))
        value_right = DCSVAL(b, nbaux-1, break_matrix_w(index_a+1,ii_b:nb), coeff_matrix_w(index_a+1,:,ii_b:nb))

        slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
        v0_fun_welfare = value_left + slope * (a - agrid(index_a))    
        RETURN
    elseif (index2.GT.0) THEN
        if (b.GE.bgrid(index2+1)) THEN
            ii_b=index1+1
            nbaux =nb-ii_b+1
            value_left  = DCSVAL(b, nbaux-1, break_matrix_w(index_a,ii_b:nb), coeff_matrix_w(index_a,:,ii_b:nb))

            IF (index2+1.EQ.nb) THEN
                value_right = v0_matrixA_w(nb,index_a+1)
            else                
                ii_b=index2+1
                nbaux =nb-ii_b+1        
                value_right = DCSVAL(b, nbaux-1, break_matrix_w(index_a+1,ii_b:nb), coeff_matrix_w(index_a+1,:,ii_b:nb))
            endif
            
            slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
            v0_fun_welfare = value_left + slope * (a - agrid(index_a))    
            RETURN
        elseif (b.LT.bgrid(index2+1)) THEN
            v0_fun_welfare = -1d6
            RETURN
        endif
    endif
elseif (index1.GT.0) THEN
    if (b.GE.bgrid(index1+1)) THEN
        if (index2.EQ.0) THEN
            
            IF (index1+1.EQ.nb) THEN
                value_left = v0_matrixA_w(nb,index_a)
            else                
                ii_b=index1+1
                nbaux =nb-ii_b+1
                value_left  = DCSVAL(b, nbaux-1, break_matrix_w(index_a,ii_b:nb), coeff_matrix_w(index_a,:,ii_b:nb))
            endif           

            ii_b=index2+1
            nbaux =nb-ii_b+1        
            value_right = DCSVAL(b, nbaux-1, break_matrix_w(index_a+1,ii_b:nb), coeff_matrix_w(index_a+1,:,ii_b:nb))

            slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
            v0_fun_welfare = value_left + slope * (a - agrid(index_a))  
            RETURN
        elseif (index2.GT.0) THEN
            if (b.GE.bgrid(index2+1)) THEN
                IF (index1+1.EQ.nb) THEN
                    value_left = v0_matrixA_w(nb,index_a)
                else                
                    ii_b=index1+1
                    nbaux =nb-ii_b+1
                    value_left  = DCSVAL(b, nbaux-1, break_matrix_w(index_a,ii_b:nb), coeff_matrix_w(index_a,:,ii_b:nb))
                endif
            
                IF (index2+1.EQ.nb) THEN
                    value_right = v0_matrixA_w(nb,index_a+1)
                else                
                    ii_b=index2+1
                    nbaux =nb-ii_b+1        
                    value_right = DCSVAL(b, nbaux-1, break_matrix_w(index_a+1,ii_b:nb), coeff_matrix_w(index_a+1,:,ii_b:nb))
                endif
                
                slope = (value_right - value_left)/ (agrid(index_a+1) - agrid(index_a))
                v0_fun_welfare = value_left + slope * (a - agrid(index_a))    
                RETURN
            elseif (b.LT.bgrid(index2+1)) THEN
                v0_fun_welfare = -1d6
                RETURN
            endif
        endif
    elseif (b.LT.bgrid(index1+1)) THEN
        v0_fun_welfare =-1d6
        RETURN
        
    endif
endif
    
end function v0_fun_welfare
!*********************************************************************************
SUBROUTINE HIST(distAaccess,distBaccess,distAaut,distBaut)

USE GLOBAL

IMPLICIT NONE

REAL(8),DIMENSION(nb2,na),INTENT(OUT)::distAaccess,distBaccess
REAL(8),DIMENSION(na),INTENT(OUT)::distAaut,distBaut
INTEGER::maxdata,i,i_a,i_b,nbA,nbB
PARAMETER(maxdata=10000)
REAL(8), DIMENSION(maxdata)::dataA,dataB
REAL(8)::sum1,sum2

nbA = nb2
nbB = nb2

OPEN(1,FILE='output//welfare//mcdistA2.txt')  
OPEN(2,FILE='output_baselineX\welfare_uncond\mcdistA501.txt')   

DO i=1,maxdata
    READ(1,*) dataA(i)
    READ(2,*) dataB(i)
ENDDO

CLOSE(1)
CLOSE(2)

distAaccess = 0d0
distBaccess = 0d0
distAaut = 0d0
distBaut = 0d0
sum1=0d0
sum2=0d0

DO i_a=1,na
    DO i_b=1,nbA
        distAaccess(i_b,i_a) = DBLE(count(dataA.EQ. i_a +(i_b-1)*na))/DBLE(maxdata)
        sum1=sum1+distAaccess(i_b,i_a) 
    PRINT*,i_a,i_b,distAaccess(i_b,i_a) 
    ENDDO
    DO i_b=1,nbB
        distBaccess(i_b,i_a) = DBLE(count(dataB.EQ. i_a +(i_b-1)*na))/DBLE(maxdata)
        sum2=sum2+distBaccess(i_b,i_a) 
    ENDDO
ENDDO


DO i_a=1,na
    distAaut(i_a) = DBLE(count(dataA.EQ. na*nbA+i_a))/DBLE(maxdata)
    distBaut(i_a) = DBLE(count(dataB.EQ. na*nbB+i_a))/DBLE(maxdata)
    sum1=sum1+distAaut(i_a) 
    sum2=sum2+distBaut(i_a)
ENDDO

IF (DABS(sum1-1d0).GT.1d-5) THEN
    PRINT*,'Probs of Alternative do not sum up to 1',sum1
ENDIF

IF (DABS(sum2-1d0).GT.1d-5) THEN
    PRINT*,'Probs of Baseline do not sum up to 1',sum2
ENDIF


END SUBROUTINE HIST
!*****************************************************************************
SUBROUTINE HIST2(distAaccess,distBaccess,distAaut,distBaut,debtyT_A)

USE GLOBAL

IMPLICIT NONE

REAL(8),DIMENSION(nb2,na),INTENT(OUT)::distAaccess,distBaccess
REAL(8),DIMENSION(na),INTENT(OUT)::distAaut,distBaut,debtyT_A
INTEGER::maxdata,i,i_a,i_b,nbA,nbB,bp1_A(na,nb2)
PARAMETER(maxdata=10000)
REAL(8), DIMENSION(maxdata)::dataA,dataB
REAL(8)::sum1,sum2,debt_sum,debt_val,debt_count

nbA = nb2
nbB = nb2

open (110, FILE='output\bgrid2.txt')
DO i=1,nb2
    READ(110,*) bgrid2(i)
ENDDO

CLOSE(110)


open (3, FILE='output\welfare\gibp_501.txt')

do i_a = 1,na
    do i_b = 1,nb2
    READ(3, '(I7)') bp1_A(i_a,i_b)
    end do
end do

CLOSE(3)

OPEN(1,FILE='output//mcdistA2.txt')  
OPEN(2,FILE='output_baselineX\welfare_uncond\mcdistA_501.txt')   

DO i=1,maxdata
    READ(1,*) dataA(i)
    READ(2,*) dataB(i)
ENDDO

CLOSE(1)
CLOSE(2)

distAaccess = 0d0
distBaccess = 0d0
distAaut = 0d0
distBaut = 0d0
sum1=0d0
sum2=0d0

DO i_a=1,na
    debt_sum=0d0
    debt_count=0d0
    DO i_b=1,nbA
        distAaccess(i_b,i_a) = DBLE(count(dataA.EQ. i_a +(i_b-1)*na))/DBLE(maxdata)
        debt_val = bgrid2(i_b)
        debt_sum = debt_sum+debt_val*DBLE(distAaccess(i_b,i_a)) 
        debt_count = debt_count+distAaccess(i_b,i_a) 
        sum1=sum1+distAaccess(i_b,i_a) 
    ENDDO
    
    debtyT_A(i_a)=debt_sum/DBLE(debt_count)

    DO i_b=1,nbB
        distBaccess(i_b,i_a) = DBLE(count(dataB.EQ. i_a +(i_b-1)*na))/DBLE(maxdata)
        sum2=sum2+distBaccess(i_b,i_a) 
    ENDDO
ENDDO


DO i_a=1,na
    distAaut(i_a) = DBLE(count(dataA.EQ. na*nbA+i_a))/DBLE(maxdata)
    distBaut(i_a) = DBLE(count(dataB.EQ. na*nbB+i_a))/DBLE(maxdata)
    sum1=sum1+distAaut(i_a) 
    sum2=sum2+distBaut(i_a)
ENDDO

IF (DABS(sum1-1d0).GT.1d-5) THEN
    PRINT*,'Probs of Alternative do not sum up to 1',sum1
ENDIF

IF (DABS(sum2-1d0).GT.1d-5) THEN
    PRINT*,'Probs of Baseline do not sum up to 1',sum2
ENDIF


END SUBROUTINE HIST2
!*****************************************************************************
